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This work addresses and solves the problem of generically 
tracking black hole event horizons in computational simula- 
tion of black hole interactions. Solutions of the hyperbolic 
cikonal equation, solved on a curved spacetime manifold con- 
taining black hole sources, are employed in development of a 
robust tracking method capable of continuously monitoring 
arbitrary changes of topology in the event horizon, as well 
as arbitrary numbers of gravitational sources. The method 
makes use of continuous families of level set viscosity solu- 
tions of the eikonal equation with identification of the black 
hole event horizon obtained by the signature feature of dis- 
continuity formation in the eikonal's solution. The method is 
employed in the analysis of the event horizon for the asym- 
metric merger in a binary black hole system. In this first 
such three dimensional analysis, we establish both qualitative 
and quantitative physics for the asymmetric collision; includ- 
ing: 1. Bounds on the topology of the throat connecting the 
holes following merger, 2. Time of merger, and 3. Continuous 
accounting for the surface of section areas of the black hole 
sources. 



I. INTRODUCTION 

The idea of a totally collapsed gravitational source, 
from which nothing — not even light — can escape, is 
an old one, dating back at least to the work of Laplace 
(and others) in the eighteenth century [1]. Since those 
original considerations, research has uncovered a wealth 
of understanding regarding the physics of black holes. In 
general, the work has followed two major routes, with 
an ever - narrowing gap between the approaches. Along 
one direction, black holes are studied as mathematical 
solutions in a given theory of gravity including Newto- 
nian theory, post-Newtonian metric theories of gravity, 
Einsteinian gravity, as well as semi-classical and (more 
recently) various attempts at full quantum theories of 
gravity. The second direction considers black holes in the 
astrophysical and astronomical contexts. Both fields of 
research have seen increasing activity over the years and 
have made startling discoveries concerning the physics of 
black holes; including the gravitational collapse theory of 
Oppenheimer and Snyder [2], proofs of uniqueness and 
stability of black holes in Einstein's general relativity [3], 
thermodynamic properties of black holes including the 
three laws of black hole mechanics [4] and mechanisms for 
black hole radiation and evaporation [5] , discovery of crit- 
ical phenomena in black hole formation [6], experimen- 
tal identification of both astrophysical black hole sources 
themselves [7] and their event horizons [8], experimen- 
tal signatures for super - massive black holes in galactic 



centers [9] , and experimental bounds on the distribution 
and spectrum of black hole sources and collisions [10]. 

Numerical relativity, which addresses computational 
solution of Einstein's equation [11], is an active partici- 
pant in both the mathematical and the astrophysical pro- 
grammes of research. With the advent of the Laser Inter- 
ferometric Gravitational Observatory and other similar 
efforts, considerable attention has focused on the generic 
binary black hole coalescence (BBHC) problem, expected 
as the strongest sources for gravitational waves [12]. As 
discussed in [11], this is a very difficult and interesting 
problem that can only successfully be addressed in the 
computational domain. Consequently, the binary black 
hole coalescence problem has become an active subject 
of numerical relativity. One particular problem in the 
computational domain is the computational definition, 
detection and tracking of the black hole event horizon 
itself [13], [14], [15], [16], [17], [18]. 

The present work completely solves the problem of nu- 
merically tracking black hole event horizons. The so- 
lution is complete in the sense that a single method 
is presented such that any one implementation of the 
method can generically detect arbitrary numbers of black 
hole event horizons undergoing arbitrarily strong gravi- 
tational interactions. For example, using a single compu- 
tational code of our method we analyse both single black 
holes and black holes undergoing merger; and no special 
modifications of our code are required to handle these 
distinct dynamics. 

The present article, describing our generic method for 
tracking black hole event horizons, is divided as follows. 
In section II the eikonal equation, the foundation of our 
method, is described in sufficient detail to be employed in 
an event horizon tracker. In particular, we focus on the 
signature behavior of a black hole event horizon in solu- 
tions of the eikonal equation. As usual in numerical work, 
there are a variety of possible implementations of the ap- 
proach. In section III, several closely related systems of 
equations are presented and one particular system is sin- 
gled out for consideration. The system chosen makes use 
of an explicit second order diffusion, or viscosity, term 
and we show in section III the relationship between so- 
lutions of the diffusive equations of motion and the con- 
tinuum cikonal equation of interest. Section IV details 
extraction of the two dimensional sections of the event 
horizon for each time level of a numerical evolution. Such 
extraction is crucial for carrying out area, mass, and spin 
calculations. Section IV also shows the accuracy of our 
implementation by considering a parameter space survey 
of single Kerr black holes. Sections V - X describe in 
detail the first three dimensional application of our (or 
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any) method to the numerical analysis of the asymmet- 
ric binary black hole coalescence problem. In particular 
we search for evidence of any nontrivial topology of the 
horizon immediately following merger. In distinction to 
the prediction of [23] we find a topologically spherical 
horizon to within the accuracy of our three dimensional 
mesh. Further, area analysis of the candidate numerical 
event horizon is carried out in conjunction with analysis 
of the black hole apparent horizons, which reveals both 
a time of merger much earlier then estimated using ap- 
parent horizons and mass energy estimates much larger 
those found using apparent horizon tracking methods. 
Our conclusions are presented in section XI. 



II. THE EIKONAL 

To begin, since the Lagrangian L = g ab x a x b of null 
geodesic motion has only kinetic terms it is equal to the 
associated Hamiltonian 1 . Legendre transformation 

H(T, X \p b ) = ^p c -LUx d ^\ (1) 



where 



Pa 



dL 



3©' 



(2) 



sets L = H. The corresponding Hamiltonian equations 
are 



dx a dH „ . 

dT Opa 



dpb 
dT 



— = -d b H = -pcPdd b g 



cd 



(3) 



(4) 



It is generic that the spacetime metric is independent 
of the affine parameter r. There is thus a first integral 
associated to geodesic motion. Making use of this prop- 
erty permits elimination of the affine parameter in favor 
of coordinate time t. To see this, it is convenient to adopt 
the ADM variables 



9 tt = ^,9 U = ^,9 ij =l ij -^- 
or a z a z 



(5) 



Here 7 y is the inverse spatial metric. Substituting the 
ADM variables directly into the Hamiltonian gives 



H = Pa g ab Pb 



(6) 



Since the Hamiltonian is not explicitly r dependent there 
is a constant of the motion 

-^Pt + + W (V - Pi = - 2 - (7) 

For uj 2 < 0, uj 2 = 0, or u 2 > the motion is said to be 
timclike, null, or spacelike. Without loss of generality, as- 
suming null geodesic motion, solution of (7) by ordinary 
algebra yields: 



Pt = (3 l pi ± ayJpiYipj. 



(8) 



With this result, the Hamiltonian can be explicitly fac- 
tored: 



H = H + H- = 



where 



H±=p t - I3 l pi ± a^pa^pj = 0. 



(9) 



(10) 



In the case of either root Hamilton's canonical equations 
become 



dt dH ± 



dr 



dpi 



\/pk-f k 



Pi 



dp t = H dH± 



dT 



dpi 
dT 



= -H-. 



dt 



dH± 



(11) 



(12) 



(13) 



(14) 



According to equation (11), r can be eliminated in favor 
of t. The system of equations becomes simply the equa- 
tions (12) and (14) with t written in place of r and the 
factor H T cancelled everywhere on the right hand side. 
Explicitly, 



dt 



\/pkl k 



(15) 



Pi 



and 



dpi 
dt 



= -di (-P Pj ± a^ifS] • (16) 



The eikonal, corresponding to the Hamiltonian of equa- 
tion (6), 



1- We use early- latin indices to denote spacetime components 
a,b,c — 0, 1, 2, 3; mid-latin indices denote spatial components 
i,j,k = 1,2,3. 



d a Sg ab d b S = 



(17) 



can be factored similarly. To do so it is a simple matter 
of making the replacements p t — ► d t S, pt — ► diS in (8) to 
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find the following symmetric hyberbolic partial differen- 
tial equation 



d t S = (3 l d t S ± a^jdiS-fidjS = H. 



(18) 



Note that a bar is introduced here to distinguish the 
Hamiltonian used here from the Hamiltonian used in (6). 
This equation is also used in the method of [13], although 
in that work the equation is further reduced to consider 
the case of a single null surface. 

The right hand side of this result is homogeneous of 
degree one in diS. The characteristic curves along which 
the level sets Y of S are propagated, are then 



x l = -f3 i ±a 



d l S 



dH 



d t S = diH(t, x>,djS) 



(19) 



(20) 



which are the null geodesies of equations (15) and (16). 
Immediately, the integral curves of the gradients of S and 
T are also the null geodesies: 



dx l (A) 
dX 



g ia Pa = g ta d a S = d l S (A, x j (A)) . (21) 



Hereafter, the bar on H is dropped with the under- 
standing that the Hamiltonian considered is that of equa- 
tion (18). This result establishes that the eikonal is tech- 
nically a Riemann invariant of the null geodesies, a fact 
that proves useful in establishing the signature of a black 
hole event horizon in solutions of the eikonal equation. 
More specifically, since the eikonal is the canonical gen- 
erator of null geodesies, it can be employed in analysis of 
black hole event horizons, which are by definition gener- 
ated by null geodesies having no future end points. To 
proceed in this manner the equation (18) is employed in 
an initial value problem and then surveyed for signature 
features of black hole event horizons. 

A black hole event horizon is generated by a congru- 
ence of outgoing — but future asymptotically nonex- 
panding — null geodesies. The scope of the surveys of 
the eikonal equation that are required to identify a black 
hole event horizon is then restricted to the space of all 
outgoing null surfaces. These surveys are greatly reduced 
by the fact that solutions of the eikonal are categorized 
by topologically equivalent solutions. To see this, note 
that 



and 



d i ip(S) = ^d i s = \(S)d i s, 



d t il>(S) = ^d t S = \(S)d t S. 



(22) 



(23) 



By homogeneity of the right hand side of (18), if S is 
a solution then tp (S) is also solution. Thus smoothly 



related initial data So — > S' — tp (So) have smoothly 
related solutions. This feature alleviates the need for 
surveying over smoothly related initial data. 

A further reduction of the scope of solution surveys 
is provided by the equivalence of ingoing and outgoing 
solutions under time reversal. Propagation of data for 
S describing an ingoing or outgoing null surface is ac- 
complished by specification of: 1. A definition of the 
direction of time, 2. A choice of a and /3 l , and 3. A 
choice of the root. With these choices specified, data is 
then uniquely partitioned into an ingoing type and an 
outgoing type with the distinction being the gradient of 
S. 

With the dynamics of any outgoing null surface (in- 
cluding the event horizon) specified by (18), and the 
scope of solution surveys categorized into topological 
classes, the task remains of identifying the event horizon 
within this restricted space of outgoing null surfaces. To 
do so, the results and approach of [13] are adopted here 
and modified to include the eikonal equation as discussed 
above. The result of this approach is a signature feature 
of black hole event horizons in the eikonal equation. 

An event horizon of a black hole is by definition a criti- 
cal outgoing surface when tracked into the future. Let V 
be a point interior to the horizon and Q be a point exte- 
rior to the horizon such that V and Q lie on characteris- 
tic curves of the eikonal 77? and jq. At arbitrarily early 
times let 77? and 7q pass arbitrarily closely to a point Ti 
that lies on the horizon Y. Since S is a Riemann invariant, 
at arbitrarily early times the jump of the eikonal at Ti be- 
comes [[S {H)}} = S- (Ti) - S+ (Ti) t*S(V)-S (Q) ^ 0. 
In the computational domain, where the resolution is fi- 
nite, this discontinuity will appear generically in finite 
time. As such, an approximation of the event horizon 
will appear numerically as the formation of a jump dis- 
continuity in the eikonal for outgoing data that is prop- 
agated into the past. This is the numerical signature of 
black hole event horizons in the eikonal's solutions. 



III. NUMERICAL METHODS 

Analysis of the continuum properties of the eikonal 
equation as a Hamilton Jacobi equation identifies three 
closely related approaches to tracking black hole event 
horizons: 

• System I: The null geodesic equations 

S {x j ,t) = S (x k ) - J dtH (t,x\ Pl ) (24) 

where the integral is evaluated along the solutions 
of 



x i = i ±a 



P' 



(25) 



-diH. 
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• System II: The eikonal 

d t S = -H(t,x j ,d j S) (26) 

• System III: The flux conservative form 

S(x\t) = <j>dx l p t (x j ,t) (27) 

dm (x\t) + d t H (t, x j , Pj ) = 0. (28) 

In each of the above three systems of equations the 
Hamiltonian is given by (18). 

In each of the above cases, the structure of the dynam- 
ical equations provide certain advantages. For example, 
in each case the symplectic structure can be employed to 
identify numerical loss of accuracy in a similar manner 
to some modern numerical schemes used in Hamiltonian 
dynamics. Further, as a flux conservative system, high 
resolution methods from computational fluid dynamics 
can be applied directly to the third system. In following 
sections of this article, the second system is considered 
since this system of equations yields an expedient imple- 
mentation that is sufficiently accurate for our purposes. 

Singular behavior on the eikonal is not specific to only 
the event horizon and instead, as described in detail by 
Arnold and Newman [19], [20], [21], the eikonal is known 
to generically break down on caustic and other sets. Spe- 
cial numerical methods arc then required to handle the 
generic singular behavior of the eikonal; we make use of 
an explicit viscosity term. In the continuum, addition 
of our form of numerical viscosity at the level of the fi- 
nite difference approximation corresponds to replacing 
the evolution of 



d t S = -H(t,x i ,d j S) 
with evolution of the equation 

dfip = e 2 V 2 V - H(t,x\dji>) 



(29) 



(30) 



where e is a small quantity which we call the viscosity and 
V 2 denotes any second order, linear derivative operator. 
There is a well defined sense in which the solutions S 
relate to the solutions ip; it is simply given (when the 
solutions S exist) by the WKB transformation 

iP(x\i) = ^ a„ (xV) e™ exp (S7e) = Aexp(S/e). (31) 



To see the explicit relationship between the solutions ip 
and the solutions S note that 



d t i> = ^{d t S + ed t logA), 



d i iP = ^(d i S + ed i logA). 



(32) 



(33) 



V 2 V = J (VS + eV log Af + ^(W 2 S + eV 2 log A) . 

(34) 

Assuming that the Hamiltonian is homogeneous of de- 
gree one in momentum, and making use of perturbation 
theory: 



H(x\Vjip) = ^ (H(x i ,V j S)+eH 1 (t, x l , dj log A) ) . 

(35) 



Substituting these results into (30) and cancelling an 
overall factor of ip/e gives at lowest, and first order: 

d t S = -H(t,x i ,d j S). (36) 

ed t log A = -eHi (i, x l , dj log A) (37) 
+ e(V5' + eVlogA) 2 + ... 



Here H± is the first order linear Hamiltonian obtained 
from perturbation theory. At zeroth order, S then satis- 
fies the eikonal equation, while the first order correction 
is the linear result of first order perturbation theory and 
expresses the evolution of A. Higher order results can 
be found similarly. Again, these results hold only where 
solutions of S exist; that is, away from discontinuities or 
other solution singularities in S and its derivatives. 

The above analysis is a slight modification of usual 
WKB expansion and expresses the solutions ip in terms 
of the solutions S. Such solutions can be inverted using 
the conventional method of series inversion. To invert 
the WKB expression, and express the solutions of the 
eikonal in terms of the solution of the parabolic equation, 
assume first that the solution of the parabolic equation 
(30) is given as ip (numerically or otherwise). Writing 



(38) 



the procedure is then analogous to that of the WKB ex- 
pansion: A is constructed as an asymptotic power series 
with coefficients depending on ip alone. The result of 
such an analysis is simply: 



S(x,t) = log ip(x,t) + e / dt 



,V 2 V(x,i) 



4> 



where the integral is evaluated along 

dx l dH (t, x j , d S k log tp) 



dt 



dpi 



(39) 



(40) 



and ip is provided independently by (numerical) solution 
of 



d t ip — e\7 2 tp — H (t, x\djip). 



(41) 
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One advantage of an explicit second order viscosity 
term is a simple procedure for reducing the error of vis- 
cosity solutions by one order of the viscosity e. The limit 
e — > corresponds to a continuous family of zeroth order 



solutions 4> t ] where 



= S + Ve 



(42) 



and rj e = O (e). Note that each solution cf> e is obtained 
from an analysis similar to that following equation (38), 
although obtained by neglecting the logarithm. 

Given two viscosity solutions (j) e and </> 2 e it is then pos- 
sible to construct a third improved viscosity solution <f>i 
that is accurate to O (e 2 ). To see this, consider the com- 
bination 

4>! = <j>, + {<t> e - <j> 2e ) . (43) 
4>i is accurate to O (e 2 ) since 

4>e ~ <j>2e =Ve~ V2e ~ Ve ~ ^Ve = ~Ve- (44) 

To make use of these improved viscosity solutions let h 
denote the resolution of the numerical mesh. Any second 
order finite difference approximation will have then 



= 4> + o(h 2 ), 



(45) 



where <j) is the continuum solution and <f> denotes its finite 
difference approximation. [We will use a hat throughout 
the paper to denote the discrete approximation to a con- 
tinuum object.] Similarly, 



Using 



gives 



S = S + O (h 2 ) . 



S = ^ + 0{e) = <f> I + 0(e 2 ) 



(46) 
(47) 



S = $ + O (h 2 ) +0{e) = fa + O {h 2 ) + O (e 2 ) . (48) 



where r = u(6, cj>) is the surface function for a given cen- 
ter c l contained within the surface V. Given a choice for 
the center, the radial function for S (T) = S = can 
then be approximated via the interpolation and bisec- 
tion. 

To establish the accuracy of our implementation, in- 
cluding the routines that accomplish extraction of the 
level sets and the accuracy of the viscosity term, we con- 
sider stationary, spinning black holes. This case is com- 
pletely described by the Kerr - Newmann family of ax- 
isymmetric solutions of Einstein's equation. It is conve- 
nient to make use of the Kerr-Schild form for the metric 
[25] since this form is used in our binary black hole evo- 
lution code as well as in our solution of the initial value 
problem for setting initial data for the evolution. Specif- 
ically, the metric in the Kerr-Schild form is 



g ab = ^ab _ 2Hl a l b . 



(49) 



Here r\ a b = diag (—1, 1, 1, 1) is Minkowski's metric, H is 
a space time scalar, and l a is an ingoing null vector with 
respect to both the Minkowski and full metric. The Kerr 
solution is the two parameter family of solutions such 
that 



H 



Mr 3 



2,2 



and 



r 



V = 



rx 



ay 

2 ' 



H + a 



ry — ax 
r 2 + a 2 ' 



(50) 

(51) 
(52) 

(53) 
(54) 



IV. LEVEL SET EXTRACTION 

Of crucial interest in the binary black hole coalescence 
problem are the areas of sections of the black hole event 
horizon. To find these sections, at any given time level 
any of the level sets, say 5 = 0, can be extracted to ob- 
tain the surface T, a two dimensional section of the cor- 
responding null surface. This problem of extraction is an 
inverse problem, since it requires that points (x,y, z) are 
found such that S (x, y, z) = 0. To accomplish this inver- 
sion, we find that an ordinary bisection method is suffi- 
cient for use with an ordinary second order interpolation 
scheme. In this method it is assumed that the surface Y 
can be expressed in spherical coordinates (8, </>, u (9, 0)), 



r 2 = \ (p 2 a 2 ) + yjl {p 2 -a?)+a?z 2 (55) 

where 

p 2 = x 2 + y 2 + z 2 . (56) 

Finally, the event horizon for the Kerr black hole is lo- 
cated on the ellipsoid r = r + = r (x,y, z) where 

r+ = M + VM -a. (57) 
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t (Units otM) 



t (Units olM) 



FIG. 1. Percent error in area of M — 2, a = event hori- 
zons in survey over viscosity parameter: e = h 2 , h 2 /2, ft 2 /4, 
h 2 /8. Here increasing £ corresponds to propagation into the 
past. 

In figures (1) -(4) we show the evolution backward in 
time of the eikonal equation (followed by extraction of the 
S = surface). Errors can arise both in the evolution 
and in the extraction of the surface. Figure (1) shows 
the percent errors in the extracted areas for a nonspin- 
ning black holes with mass M = 2,a = 0ina survey over 
the viscosity parameter. Figure (2) shows a similar study 
but for the L2 norm of the truncation error in the func- 
tion r + . Note that the function r + is defined for every 
point of the horizon in the continuum r+ = r+ (x,y, z). 
In the computational domain r + then takes a discrete 
form f + = f + k), where the integers k span the 
numerical mesh. The truncation error e r+ is then 

e r+ (i, j, k) = r + (x, y, z) - f+ (i, j, k) (58) 

where is it understood that both r+ and r + are evaluated 
at the same point. Bith Figures(l) and (2)show that 
viscosity parameter /i 2 /8 adequately captures the horizon 
location. While not perfect, it will suffice for the short 
term horizon tracking reported here. 



FIG. 2. L2 norm in truncation error of r+ for M = 2, 
a — event horizons in survey over viscosity parameter: e = 
h 2 , h 2 /2, h 2 /4, h 2 /8. Here increasing t corresponds to prop- 
agation into the past. 

However these results also suggest that in vanishing 
viscosity the percentage error in the calculated area is 
reduced toward a bias. This bias is partly associated 
to the finite resolution of the computational mesh and 
partly to accuracy of the extraction routine. These fig- 
ures were generated using a three dimensional computa- 
tional domain of N 3 points with N = 121. The outer 
boundaries are located at [— 15M, +15M] in the x,y,z 
directions. The resolution of the finite difference mesh 
for these results is then h = M/A. Also, a Courant - 
Friedrichs - Lewy factor of A = 1/4 with an iterated 
Crank Nicholson scheme [22] was used as the finite differ- 
ence approximation of the evolution of the eikonal equa- 
tion (18). Neumann boundary conditions diS = on 
the outer boundary arc found to be generically sufficient 
conditions for stability of the method. The philosophy 
here is that the primary interest is deep within the bulk 
of the computational domain where the event horizon of 
the black hole is located. The outer boundary is then 
treated only to the degree that the evolution of the in- 
terior region remains stable. Further, the interior of the 
black hole is excised from the computational domain in a 
sphere of radius r = r + 2dx, where ro denotes the radius 
of the Kerr - Newman ring curvature singularity. Finally, 
the discrete surface T was constructed using m 2 points 
with m — 100. At such a resolution of the extracted sur- 
face any errors in the area must be attributed to all of: 
the viscosity parameter, the extraction routine, and the 
resolution of the underlying three dimensional grid. 



6 




10 15 20 

t (Units olM) 



FIG. 3. Percent error in area for three dimensional level 
set solutions a/M =, 0, 1/4, 1/2, 3/4. Here increasing t 
corresponds to propagation into the past. 

Figures (3) and (4) show the percent errors in the ex- 
tracted areas as well as the L2 Norm of the truncation 
error in the function r + in a survey over the angular mo- 
mentum parameter a. Here the viscosity is e — h 2 . These 
figures both show the evolution from a single null sphere 
that is completely exterior to the horizon and propagated 
into the past. Since the event horizon of a spinning black 
hole is elliptical in its geometry, the spherical data we 
have chosen for t± = corresponds to a percent error in 
the area and r+ that varies with the spin parameter a/M 
at ij_ = 0, explaining why the curves do not intersect at 
t± = 0. (Where appropriate we append a j to t, thus: 
tj, indicating evolution into the past; we also sometimes 
use to emphasize that we mean the forward evolving, 
usual, time t. Thus t± = corresponds to the late time 
at which we begin to integrate into the past.) 

However, the errors should converge to a constant, 
which is evident in each of the curves with a < 0.75 in 
figures (3) and (4). For a = 0.75 the curve does not con- 
verge, and we expect that different choice of initial data 
will exhibit convergence. According to these results, the 
viscosity level set method does indeed accurately and ro- 
bustly detect the distorted outermost event horizons of 
spinning black holes at least when a < 0.75. Note that 
we study both the accuracy of r + and the accuracy of 
the calculated areas since the calculations seperately and 
together establish the accuracy of our area calculation 
and of our detection of the event horizon. 




5 10 15 20 25 

t (units of M) 



FIG. 4. L2 norm of truncation error in r+ for three di- 
mensional level set solutions a/M — 0, 1/4, 1/2, 3/4. Here 
increasing t corresponds to propagation into the past. 



V. ASYMMETRIC BINARY BLACK HOLE 
COALESCENCE 

Analysis of the event horizon for the binary black hole 
coalescence problem in the case of head on collision has 
been considered in detail in [16]. The problem of the 
event horizon for asymmetric, that is off axis, collision 
has only been considered analytically [23] and prior to 
this work no results for numerically generated sources 
have been analyzed. Numerical evolution and analysis of 
an asymmetric binary black hole system was studied in 
[24], but at that time the question of the event horizon 
was not considered. 

In this section the results of the previous sections are 
applied to the first completely numerical analysis of the 
event horizon for the case of asymmetric collision. 

To begin, consider two black holes of mass M = 1 with 
aligned spins in the positive z direction of a/M = 1/2. 
The computational domain is a grid of A 3 points with 
A = 121. The outer boundary is located at ±15M and 
the holes are initially positioned at (x, y, z) = (+6, +2, 0), 
and (x,y,z) = (—6,-2,0). This computational domain 
is identical to the mesh used in the previous section to 
analyze the percent error in area calculations of surfaces 
extracted from the level set method. The percent error 
in the calculation of the area of sections of the horizon 
should have a magnitude of about 4 — 5%. Further, as an 
order of magnitude estimate, in a flat spatial geometry; 
e.g., in a Newtonian spacetime, the initial separation s of 
the black hole centers would be s = \J\2 2 + 4 2 w 12.64. 
This would seem to be an ample initial separation to 
guarantee that the initial data corresponds to two dis- 



7 



tinct black holes. However, in the case that the holes are 
nonspinning, each will have a spherical event horizon of 
radius r = 2M. Assuming only marginal distortion of 
the nonspinning event horizons due to spin effects, which 
is an approximation that is justified by the properties of 
spin a/M — 1/2 black holes, the nearest separation be- 
tween the two sections of the black hole event horizon 
is then approximately s min = 12.64M - 4M w 8.64M. 
Again, this approximation assumes a flat underlying ge- 
ometry; and so can be considered only as an order of 
magnitude estimate. These initial data then appear to 
correspond to a separation of approximately two nonspin- 
ning black hole diameters between the surfaces of section 
of the black holes. 

The holes are boosted along the x direction with speeds 
of ±c/2. This boost lengthens the nearest separation 
Smin of the holes due to Lorentz contraction of the hori- 
zons. [In this coordinate system the horizons undergo con- 
traction in the direction of motion. For a single hole, the 
area of the horizon does not change under this boost.] 
The nearest coordinate separation s TO j„ between the holes 
is then expected to lie in the range 8.64 < s min < 12.64. 

The numerical evolution of this collision process is 
carried out for approximately 10 M of run time with a 
Courant factor of A = dt/dx = 1/4. The code is the 
Texas black hole evolution code, a derivative of the Agave 
code [24]. Apparent horizon finders [25] locate two dis- 
tinct apparent horizons of area A w 50M for the initial 
data and continue to do so until t = 8M, when only a sin- 
gle apparent horizon of area A w 200M can be located. 
This single apparent horizon persists until approximately 
10M, beyond which instability effects, stemming from 
the outer boundary and the excision boundary, swamp 
the solution. 



VI. DATA FOR THE EIKONAL 

While the run length of this asymmetric collision data 
is relatively small in units of the black hole masses, the 
problem of detecting the associated black hole event hori- 
zon, or horizons, is not a small computational problem. 
For example, to analyze this data requires analysis of 
the lapse a, the shift vector and the three met- 
ric 7 1 - 7 . By symmetry of the metric, 7 y = 7^, there 
are only 6 independent components. Tracking the event 
horizon of the associated data then requires analysis of 
10 grid functions, where each grid function consists of 
O (Nt x N 3 ) x 8B = O (160 x 121 3 ) x 8B w 2.26GB 
of data. That is, tracking the associated event horizon 
requires analysis of approximately 20GB of data. 

A more serious difficulty associated to this data set is 
the relaxation time, t « AM, that is typically required for 
outgoing data to converge onto the event horizon when 
followed into the past. Assuming that the collision time 
is (as suggested by the apparent horizon solvers) near 
t = 8M, perturbation theory implies that the resulting 



horizon should undergo quasi normal ringing for another 
t w 20M from that time. That is, at the time level 
t = 10M, where the event horizon tracker will begin 
tracking into the past, it is expected that the event hori- 
zon remains highly distorted and far from its stationary 
regime. The problem is to determine good data for the 
cikonal at the time level t = 10 M, which can then be 
propagated into the past. In contrast, in analysis of the 
event horizon for the case of head on collision, researchers 
made use of approximately 100M of data [16], and as- 
sumed that the final state at t = 100M was a stationary 
or quasi stationary black hole. In such a circumstance 
using spherically symmetric data at t = 100M for an 
event horizon tracking method should prove sufficient. 
The difficulty with a data set only of length t = 10M is 
that using spherically symmetric data for the eikonal at 
the time level t = 10M does not accurately approximate 
sections of the event horizon at that time. 

To accommodate data of this 10M length and to gen- 
erate sufficient data for the eikonal at the time level 
t = 10M, a method for finding a candidate section of the 
event horizon at t = 10 M is proposed here and applied 
to the problem. This method makes use of the analytic 
properties of apparent and event horizons. 

To motivate the method, note that if the event horizon 
were stationary at t = 10M then the expansion of the 
surface would be vanishing there: 

1 HA 

0(i=lOM) = -— =0. (59) 

In such a circumstance the apparent horizon could be 
used as initial data for the eikonal, which could then 
be propagated into the past. However, according to the 
second law of black hole mechanics, in the case that the 
horizon is nonstationary, which is the situation expected 
for this asymmetric problem, at t = 10M the horizon will 
satisfy 

6{t= 10M) > 0, (60) 

in the forward time direction. In the backwards time 
direction these dynamics correspond to 6 < at t = 
10M. As such, the Taylor series of any compact null 
surface, including the critical surface of the event horizon, 
is at least of the form 

A(t + dt) = A(t) + dtA8 + dt 2 ^ + ... (61) 

and probably cannot be truncated to below 

A (t + dt) w A (t) + dtA6 + ... (62) 

The apparent horizon is then a poor estimate for the 
the event horizon at t = 10M. The objective is to use 
this behavior (62) in combination with the property that 
outgoing null data followed into the past converges to 
the event horizon. The hope is to establish a better ap- 
proximation for the structure of the section of the event 
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horizon at t = 10M, which can then be used in an event 
horizon tracking method. 

To proceed, consider three compact null surfaces of 
outgoing data: One completely interior to the horizon 
r 4 , one completely exterior to the horizon r e , and one 
that is a surface of section of the event horizon T h . Let 
i| — * oo denote propagation into the past. By the prop- 
erty that outgoing data numerically converges to the 
horizon when propagated into the past, limt^ooF 1 = 
linit^oor 6 = lim ti ^ oo r' 1 . Thus in most cases, if the 
surface T l (ti + dti) is pulled back to the time level tj_ 
and compared to the surface F* (i j ) it will be completely 
exterior to T l (ijj. Similarly, if the surface T e (ti + dti) is 
pulled back to the time level £j_ it will be completely inte- 
rior to the surface T e (tjj. That is, by iteratively pulling 
the surfaces back to a single time level, spherically sym- 
metric data will approach the numerical event horizon. 
Exceptions to this general behavior stem from the pres- 
ence of caustics in the spacetime, where null surfaces in- 
tersect, and in neighborhoods of the event horizon. In 
particular, as discussed above, the horizon of this asym- 
metric collision data is expected to satisfy 8 > 0. If the 
surface T h (ti + dt±) is then pulled back to t± it will be 
interior to the surface T h Let Y\ (t±) be a pertur- 
bation of T h (ti) that is arbitrarily close to the horizon 
T h (ti) but completely interior to the surface T h (ti). Due 
to round off and truncation error of the finite difference 
approximation and the fact that lim^-,,^ T\ = T h , if the 

numerical finite difference approximation T l e (t^ + dti), 
which approximates the continuum section (ij + dti), 
is pulled back to the time level it will generically in- 
tersect L* (ti) and contain neighborhoods that are both 
interior and exterior to the finite difference approxima- 
tion r* (ti). Similar behavior will hold for numerical data 
Tl (ti) defined to be a perturbation of T h (ti) that is ar- 
bitrarily close to the horizon T (tjj , and completely 
exterior to T h (ti). 

According to this behavior, data for the eikonal at 
the time level t — 10M can be constructed by itera- 
tion over that time slice. This approach is similar to 
treating the event horizon as if it were an apparent hori- 
zon, although modified to account for this nonstation- 
ary regime. Beginning with spherically symmetric initial 
data that is well exterior to the horizon ff (i^ = 10M) 
the data is updated, pulled back to the original time 
level and reset as follows: T% (tj.) = (tj. + dti). The 
step r^ +1 (ti) — Lfj (ti + dti) i s tnen repeated for sev- 
eral hundred iterations, which corresponds to O (10) e 
folding times. Note that in the case that the spacetime 
is stationary the surface T will converge to the appar- 
ent horizon according to this method. Since the appar- 
ent horizon of a stationary spacetime coincides with the 
event horizon, this method will generically find the event 
horizon in the case of stationary spacetimes. However, in 
the nonstationary regime, which is the case for the asym- 
metric collision problem, the result of this procedure is 



at best an improved initial guess for an event horizon 
tracker. Further information, such as area analysis or 
study of the apparent horizon, is required to argue that 
the resulting surface is a candidate section of the event 
horizon. 




FIG. 5. Tanh data (see Eq(63) for the eikonal in asymmet- 
ric binary black hole coalescence. The black contour is the 
estimate of the event horizon section, while the wire mesh is 
the apparent horizon. These data are set at the end of the 
computational evolution (t = 10M), and will be evolved into 
the past. 

We show in figure (5) the result of such a process, ap- 
plying w 200 iterations over the time slice t — 10M. Fig- 
ures (5) through (11) show the eikonal function in the 
z = plane. The location of the determined guess for 
the surface T (which we will take to begin our evolution 
of the horizon into the past) is encoded into the eikonal 
by the color map. With this T, data for the eikonal can 
be written in the form: 

5(0,0!*) = 1 + tanh^^^ (63) 

In equation (63) the first argument of the eikonal is 
ti = 0; ti will increase into the past. Also, r c denotes 
the data T and c controls its steepness. Typically we 
take a transition width c on the order of a computational 
zone. Note that this surface is not considered to be a 
true section of the event horizon, but instead is a good 
initial guess or candidate section, which after a few M 
will evolve into better approximation of a true section of 
the event horizon. By way of comparison figure (5) also 
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shows at t = 10M the final apparent horizon as a white 
wire mesh. Note that both the apparent horizon and 
these data for the eikonal are highly distorted from the 
stationary case. Figure (6) shows the resulting eikonal 
function S{x,y 1 z) after 2M of evolution into the past. 
Note that the surface T is not qualitatively changed dur- 
ing the evolution. 




FIG. 6. Data of figure (5) (asymmetric binary black hole 
coalescence) evolved backward from t — 10M to t — 8M. 



VII. SURFACE EXTRACTION AND APPARENT 
HORIZONS 

Figure (7) shows several frames of the evolution of the 
eikonal using a viscosity solution of e = h 2 (not of the 
improved viscosity form). This figure display the value 
of the eikonal function on the z = plane. Figure (7) 
shows the eikonal data as an elevation map and also via 
the color map of (5). Note in those figures that null 
surfaces interior to the event horizon undergo a change in 
topology and this topological transition is continuously 
monitored by the viscosity solutions of the eikonal. In 
figure (7) t — 0.562M is shown in the upper left-hand 
corner, t = 1.5M is shown in the upper right-hand corner, 
t = 2.5M is shown in the lower left-hand corner and 
t = 5.0M appears in the lower right-hand corner. 



Asymmetric BBH Collision: 
z = const cuts of eikonal 




FIG. 7. Change of topology in eikonal for asymmetric bi- 
nary black hole coalescence, shown as an elevation map. 

The figures (8) - (11) continue the sequence of figures 
(5) - (6), and show, for several values of if, the value 
of the eikonal in the z — plane; the location of the 
apparent horizon in the 3 dimensions (the white wire 
frame); and in the black wire frame, locations of sections 
of a candidate event horizon T c (t±) that is generated by 
evolution of the eikonal equation from data constructed 
using the method described in the previous section. In 
this context, the surfaces r c (ij) are extracted from the 
eikonal data using the technique described in section IV. 
Note that T c completely contains the apparent horizons 
throughout their evolution. This is a fundamental con- 
dition that any numerically constructed black hole event 
horizon must satisfy. To determine how these results de- 
pend on the initial data T c (t± — 0), choosing initial data 

fa (ti = 0) of the form 

u 5 {9, <j>) = u c (9, 4>) - 5, (64) 

permits survey about the data u c (9, 0), where u c (9, 4>) 
corresponds to the data T c (t[ — 0). Studies with 6 = 
M/2, M, 2M establish that the level sets r,5 penetrate 
both apparent horizons for any S > M/2. These results 
suggest that the true event horizon is contained in a do- 
main parameterized by < 6 < M/2. 
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FIG. 8. Asymmetric binary black hole coalescence: 
if = 8M. This is the same as figure (6) but the black wire 
mesh shows the estimated location of the event horizon. 




FIG. 9. Asymmetric binary black hole coalescence: 
if = 5M. Note that while the apparent horizons (the white 
wire- frame "spheres") are still well separate at t = 5M, the 
event horizon (black wire frame "peanut" already has one 
component only. 




FIG. 10. Asymmetric binary black hole coalescence: 
t T = 3M. 




FIG. 11. Asymmetric binary black hole coalescence: 
tf = 1M. Careful inspection of the event horizon (black wire- 
frame) suggests two separated components. 

Figures (12) and (13) show two views of the extracted 
level set S = at t = 2M. Note that this surface is highly 
distorted and shows the event horizon just after merger. 
Also shown in these figures are the apparent horizons for 
the two black holes. In these figures each color of the 
color map denotes a level set T of the eikonal. As such, 
each color represents a null surface. From these figures it 
is apparent that at this time there appears one innermost 
null surface that completely contains both apparent hori- 
zons. Thus, the results of the viscosity solutions suggest 
a merger time much closer to 2M then the 8M found 
with the apparent horizon trackers. [Note that in figures 
(12) and (13) the event horizon is shown as a white wire 
mesh, while the apparent horizons are shown as black 
wire meshes. This is opposite to the color scheme used 
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in figures 5, 6, 8 - 11, which was an independent study 
of the evolution as opposed to the study of the throat 
geometry considered here.] 




FIG. 12. Level set extraction for asymmetric binary black 
hole coalescence: I. 




FIG. 13. Level set extraction for asymmetric binary black 
hole coalescence: II. 



VIII. CHANGE OF TOPOLOGY 

As shown in figure (7) the viscosity solutions of the 
eikonal equation do continuously monitor a change in 



topology. In the case of asymmetric binary black hole 
coalescence it is conjectured that the level set T, which 
corresponds to a section of the event horizon, must take 
a higher genus topology at merger. To investigate this 
possibility, figure (14) shows the level set T viewed along 
the axis joining the centers of the apparent horizons. In 
that figure it is apparent that the throat function of the 
topological transition is elliptical in geometry. Studies in- 
dicate that this elliptical throat function persists for all 
null surfaces (i.e. those slightly inside or slightly outside 
our candidate event horizon) undergoing the topological 
transition. Further, for all of our computed transitions of 
the null surfaces, no higher genus topology is exhibited; 
instead, the elliptical geometry of the throat persists to 
the transition. These results suggest that if there is a 
non trivial topology in the sections of the event horizon 
as a consequence of the asymmetry of the merger, then 
that topology change is bounded to occur when the mi- 
nor axis of the ellipse is within one of our computational 
zones, or h = Af/4. 




FIG. 14. Throat function for asymmetric binary black hole 
coalescence: if = 1.562A/ . 
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FIG. 15. Area versus time for asymmetric binary black hole 
coalescence. The horizontal scale is t^, i.e. time measured into 
the past. The curves are (bottom to top) for e = 2ft 2 , e — h 2 , 
and the improved viscosity solution. 

X. BLACK HOLE AREAS 

Figure (15) shows an area versus time plot for this 
asymmetric collision. The curve with the lowest area 
is the result of a viscosity solution with e = 2h 2 . The 
curve with the second lowest area is the result of a vis- 
cosity solution with e = h 2 , while the topmost curve 
is an improved viscosity solution composed of the two 
higher viscosity solutions. According to these results it 
is immediately apparent that the viscosity solutions of 
higher viscosity show a merger time that is later in t-\ 
(and therefore prior in ij) then the merger time found 
with solutions constructed in the limit of vanishing vis- 
cosity. These results then indicate that the error (or bias) 
in the merger time of the viscosity solutions is directly 
related to the magnitude of the viscosity. More precisely, 
for the continuum merger time of t* and an approximate 
merger time of t*^, constructed using a viscosity solution 
of viscosity parameter e, the function /(e) = t* — f| 
is increasing in e. (Here all surface areas are calculated 
with m 2 points where m = 100.) 



Figure (16) shows several area versus time curves for 
initial data of the form 

u s {6,<j>)=u c (e,(t>)-5 (65) 

where the data u c (8, </>) is that obtained using the 
method of section VII. From top to bottom, the curves 
show areas for S = 0, M/4, M/2 and with a viscosity pa- 
rameter of e = h 2 . Recall that studies of the apparent 
horizons found that the true event horizon is contained 
in the domain < 5 < M/2. This survey over 5 is con- 
ducted in search for the convergence signature associated 
to event horizons. Due to the time scale of this data, the 
time scale of the dynamics, and the relaxation time scale 
of the event horizon tracking method, the signature is not 
clearly apparent. However, this study of the area curves 
does show convergence of the areas, which is expected for 
null data approaching the horizon when propagated into 
the past. The curve with 5 = M/2 shows behavior sug- 
gestive of an event horizon since the S ^ M/2 curves all 
approach that of 5 = M/2. The 6 = M/2 curve is there- 
fore considered the best candidate for the numerical event 
horizon of this study. Note that the sections of this event 
horizon completely contain the correct apparent horizons 
for all ti < 9M. Further, these S = M/2 data show a 
bifurcation time at t± s=s 8.3M, which corresponds to a 
merger time of about t-\ ~ 1.7 'M. This bifurcation time is 
detected by an algorithm that searches for any points 6, <f> 
of the surface such that u (0, (f>) < h. In the circumstance 
that a point is found such that u {9, (f>) < h, the bifurca- 
tion is expected to occur in a few more dt = M/16 in the 
ti direction. A few time levels prior to this bifurcation 
time in t[ (i.e. just after the bifurcation in t), the single 
merged horizon component has an area (the area of the 
level set) computed to be A — 148.9. Based on analysis 
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of area computations for exact solutions (figure(3), we 
anticipate an error in the horizon area of several percent. 
We conservatively assign an 8% error to the areas. Just 
prior (in t) to merger the area of each black hole is then 
A = 74.5 ± 6.QM 2 . These individual areas correspond to 
a Schwarzschild mass of about M — 1.48 ± 0.12. This 
result is a substantially larger mass for each hole then 
determined by the apparent horizon finders at the time 



as more accurate computational simulations of black hole 
interactions become available. 



level t 



0. Interestingly, studies have found that appar- 



ent horizons separated by about 10M show an increase 
in their mass due to the effects of binding energy [26] . In- 
dividual masses of about M ss 1.36 are then only a slight 
departure from studies that account for the binding en- 
ergy of the holes. Further, at the time of merger the 
holes have undergone 1.7 'M of evolution, during which 
the holes could accrete any surrounding gravitational ra- 
diation present in the initial data. The presence of such 
radiation would lead to larger masses then those found 
using apparent horizon finders at if — 0. However, it is 
important to note that due to the viscosity in the solu- 
tion the result M w 1.36 can only properly be considered 
as a lower bound on the calculated masses. The most 
significant contribution to any error in this result must 
stem from the relatively small time scale of this asym- 
metric collision data and coupling of that time scale to 
the e— folding time scale of this event horizon detection 
method. 



XI. CONCLUSIONS 

In this work we have demonstrated a relatively simple 
yet robust and (most importantly) generic solution to the 
problem of numerically tracking black hole event hori- 
zons. An implementation of our method made use of an 
explicit second order diffusion term to regulate the solu- 
tion singularities associated to caustics. As demonstrated 
by analysis of analytic sources, this term does introduce 
numerical error although we demonstrate our control over 
these effects and the resulting accuracy. But, the use of a 
second order diffusion term is not required by our method 
per se; and a variety of other approaches can be em- 
ployed. Examples of other methods for controlling break- 
down of the numerical solution include those classes of 
high resolution shock capturing numerical schemes that 
are used extensively in computational fluid dynamics for 
hyperbolic problems similar to the eikonal equation. 

The application of our new method for event horizon 
tracking method considered the asymmetric binary black 
hole coalescence problem, including a detailed analysis of 
areas of the surfaces of sections, the collision time, associ- 
ated apparent horizons, and the topology of the horizon. 
Due to the relatively short time scale of the collision data, 
our method was unable to demonstrate the signature of 
the black hole event horizon. We believe that this prob- 
lem is due to the data itself and not due to our method. 
We anticipate much more accurate and convincing results 
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